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Numerical studies of the plane symmetric, vacuum Gowdy universe on T 3 x R yield strong support 
for the conjectured asymptotically velocity term dominated (AVTD) behavior of its evolution toward 
the singularity except, perhaps, at isolated spatial points. A generic solution is characterized by 
spiky features and apparent "discontinuities" in the wave amplitudes. It is shown that the nonlinear 
terms in the wave equations drive the system generically to the "small velocity" AVTD regime and 
that the spiky features are caused by the absence of these terms at isolated spatial points. 
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' Spatially inhomogeneous cosmologies constitute almost a "terra incognita" for general relativity. In particular, the 
nature of the singularity in generic cosmologies is unknown although Belinskii, Khalatnikov, and Lifshitz (BKL) have 
conjectured that it is locally Mixmaster-like Spatially homogeneous universes' approach to the singularity can be 
either asymptotically velocity term dominated (AVTD) [^]|| or Mixmaster-like ||],|J . In general, AVTD singularities 
occur when the influence of spatial derivatives can be neglected. This means that the solution to the full Einstein 
equations eventually comes arbitrarily close to the solution to the truncated equations one obtains by eliminating 
the terms containing spatial derivatives In spatially homogeneous models, spatial derivatives in the metric create 
the spatial scalar curvature which appears as a potential in minisuperspace |4|,|| . An AVTD homogeneous cosmology 
eventually evolves toward the singularity as a Kasner solution. (The Kasner universe is an exact solution of the vacuum 
Einstein equations characterized by different expansion rates in different directions with flat spacelike hypersurfaces 
Mixmaster dynamics describes the evolution (of non-flat space-like hypersurfaces) toward the singularity as an 
infinite sequence of (approximately) Kasner epochs. The Kasner indices (which determine the anisotropic expansion 
[ rates) change whenever the influence of the spatial scalar curvature reappears. Its recurring influence means that the 
Mixmaster singularity is not AVTD. BKL claimed to prove that the generic singularity in Einstein's equations and, 
in particular, in spatially inhomogeneous cosmologies is Mixmaster-like. While this result remains controversial [Q, 



I . it provides a conjecture which can be tested by evolving spatially inhomogeneous collapsing universes numerically. 

As part of such a study of this issue p|-|lO|], we have considered the Gowdy universes on T 3 x R as a test case. 
These models are one of the classes discovered by Gowdy Jl^-|l^] as a reinterpretation of Einstein- Rosen waves juj . 
While their plane symmetry precludes local Mixmaster behavior, their nonlinear gravitational wave interactions lead 
k>( \ to an interesting phenomenology which forms the subject of this paper. The Gowdy singularity has been conjectured 
j_j ■ to be asymptotically velocity term dominated (AVTD) Jjl| and has been shown rigorously to be so for the polarized 
case ||. In addition, it has been shown that for the polarized case, cosmic censorship holds: There exists a natural 
foliation labeled by < r < oo such that regular initial data evolve without singularity for all t < oo where (in 
the collapse direction) a curvature singularity occurs at r = oo |lq ]. Here we shall use numerical studies to provide 
support for the extension of the polarized Gowdy results to generic Gowdy models. 

Gowdy cosmologies represent the simplest spatially inhomogeneous cosmological solutions to Einstein's equations. 
The vacuum case may be interpreted as the two polarizations of gravitational waves propagating in an inhomogeneous 
background spacetime ^,0- The model is plane symmetric with the propagation direction of the waves orthogonal 
to the symmetry plane. Near the singularity, the metric may be rewritten as locally Kasner with spatially dependent 
Kasner indices p2[ . (This is a non-rigorous statement of AVTD behavior.). Matter fields compatible with the Gowdy 
metric have been discussed by Carmeli, Charach, and Malin (l^jl^], and Ryan p9|. Here we shall consider only 
vacuum solutions. 

The main phenomena of interest in the Gowdy universe are the growth of small scale spatial structure and the 
approach to the AVTD regime. We shall demonstrate that both result from the nonlinear terms in the wave equation 
for P, the + polarization of the gravitational waves. These act as space and time dependent potentials whose effect is 



*e-mail: berger@oakland.edu, garfinkl@oakland.edu 



1 



to drive the system to velocity dominance. As the AVTD regime is reached, the potentials approach zero. Their spatial 
dependence means that wave amplitudes grow at some spatial points and decrease at others, eventually producing 
small scale spatial dependence in the waveforms. Non-generic behavior occurs at points where the coefficients of the 
potentials (separately) vanish identically. We note that this analysis in terms of potentials explains the observed 
episodic evolution of the system to the AVTD regime and the fact that the conditions for AVTD are reached at some 
spatial points before others. This should be considered together with the asymptotic expansion developed by Grubisic 
and Moncrief (GM) jL5| which displays an exponential decay to the AVTD limit, applicable in the last episode. 

The primary advantage of the Gowdy models as a test case for the study of inhomogeneous cosmologies is the 
existence of variables in which the dynamical equations for the wave amplitudes decouple from the constraints. Initial 
data for the wave equations may be freely specified (restricted only in that the total momentum along the direction of 
propagation must vanish) while the constraints appear as prescriptions for the construction of the background metric 
from the known solution to the dynamical equations. Thus the two principal difficulties of numerical relativity — solving 
the initial value problem and preserving the constraints during the evolution — become trivial. This decoupling between 
dynamical and constraint equations disappears in more complicated inhomogeneous models. However, studies of U(l) 
symmetric cosmologies (with only one Killing field) indicate that many features of the Gowdy phenomenology such 
as spiky small scale spatial structure persist in these more complicated models J8|,^0],[w| and, presumably, would be 
important in the dynamics of generic cosmologies. 

In Section 2, the Gowdy model on T 3 x R will be reviewed. Section 3 will contain a brief description of the 
numerical method. In Section 4, the phenomenology will be described. Section 5 will contain a unified explanation 
for the phenomenology. Discussion will be given in Section 6. 



II. THE MODEL 

The Gowdy model on T 3 x R is described by the metric |ll|J| 

ds 2 = e- x / 2 e T / 2 (-e- 2T dr 2 +dd 2 ) 

+e- T [e p da 2 + 2e p Q da d5 + {e p Q 2 + e~ p ) d8 2 ] (1) 

where A, P, Q are functions of 9, r. The sign of A differs from that in || as required for consistency with the form 
given there for the equations for A. The arguments in || did not depend on this error. We impose T 3 spatial topology 
by requiring < 9, a, S < 2ir and the metric functions to be periodic in 9. If we assume P and Q to be small, we find 
them to be respectively the amplitudes of the + and x polarizations of the gravitational waves with A describing the 
background in which they propagate. The time variable r measures the area in the symmetry plane with r = oo a 
curvature singularity. Einstein's equations split into two groups. The first is nonlinear ly coupled wave equations for 
P and Q (where , = d/da): 

P, TT - e- 2T P, ee -e 2P (Q 2 - e~ 2t Q 2 6 ) = 0, (2) 

Qitt - e- 2T Q, ee + 2 (P, r Q, T - e- 2T P,g Q, 9 ) = 0. (3) 

The second contains the Hamiltonian and ^-momentum constraints which can be expressed as first order equations 
for A in terms of P and Q: 

A, r - [P 2 + e- 2T P 2 + e 2P (Q 2 + e~ 2T Q 2 )] = 0, (4) 



\,e-2{P,eP,T + e 2P Q, e Q lT ) = 0. (5) 

This split into dynamical and constraint equations removes two of the most problematical areas of numerical relativity 
from this model: (1) The normally difficult initial value problem becomes trivial since P, Q and their first time 
derivatives may be specified arbitrarily. The only restriction, that the total 9 momentum in the waves vanishes, is a 
consequence of the requirement that metric variables be periodic in 9. (2) The constraints, while guaranteed to be 
preserved in an analytic evolution by the Bianchi identities, are not automatically preserved in a numerical evolution 
with Einstein's equations in differenced form. However, in the Gowdy model, the constraints are trivial since A may 
be constructed from the numerically determined P and Q. For the special case of the polarized Gowdy model (Q = 0), 
P satisfies a linear wave equation whose exact solution is well-known ffj]. For this case, it has been proven that the 
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approach to the singularity is AVTD ||. This has also been conjectured to be true for generic Gowdy models (JT5) 
We shall show in Section 5 how a numerical study of this model provides strong support for this conjecture. 
The wave equations (||) and (||) can be obtained by variation of the Hamiltonian |l3| 



2tt 



H= l -fd6[ir 



P + e tt q \ 





2-7T 

d9 [e~ 2r (P 2 + e 2P Q 2 )] = H K + H v . (6) 
o 



The variation of the kinetic energy term, Hk, yields equations of motion for the AVTD solution which arises when 
spatial derivatives are neglected. These equations are exactly solvable with the solution 

P = Pq + In [coshwr + cos tjj smhvr] — ► vt as t — + oo, 

e~ p ° s'mip tanhwr 

Q = Qo + t~ ; n — \ ~* as T ~~* °°> 

1 + cosip tanhwr 

tanhi>T + cosi/' 

7Tp = V ; — > V aS T — > OO, 

1 + cos?/> tanhwr 

7Tq = e p ° w sin -0 = 7Tq (7) 

given in terms of four functions of 9: ip, v > 0, P , and Q - Here Q x = Q + e~ p ° sin^»/(l + cos-0). The large r 
limits given above are for cos?/' ^ —1. The special case cosip = — 1 will be discussed later. While (0) may be solved 
for all these functions, it is convenient to note that pl| 



L+TT^e-^ (8) 

is the geodesic velocity in the target space (with metric ds 2 = dP 2 + e 2P dQ 2 ) of the wavemap Jl^,|l^,^]. If, as the 
singularity is approached, the behavior is AVTD, the true solution will approach the AVTD one. The exponential 
prefactor e~ 2r in Hv of (^|) makes plausible the conjectured AVTD singularity. However, P — > vt (for v > 0) as 
t — * oo. If v > 1, the term Va = e~ 2T e 2P (2,g in i?y can grow rather than decay as t — ► oo. In fact, if one assumes 
the AVTD behavior as r — > oo in (0) , the wave equations become |u| 

P, rr = _ e -2[l-]r [g /ja + e -2r v , V + ^-2.. 

TQ,r = e - 2 ^ T [Q'4+2«Vo T] (9) 

where ' = G?/ci#. These equations are inconsistent if v > 1 unless QJ^ = = Q'^, which does not occur in a generic 
(unpolarizcd) Gowdy model. However, in generic Gowdy models, one also expects isolated points (a set of measure 
zero) with Q' x = but Q'^ ^ 0. In Section 5, we shall see that v > 1 can persist at such isolated values of 9 where 
the AVTD conditions may not be satisfied. GM thus conjecture that the AVTD limit as t — * oo requires < v < 1 
everywhere except, perhaps, at a set of measure zero (isolated values of 9) We note that the polarized Gowdy 
model (Q = ttq = 0) has the AVTD solution P = Pq + err, irp — a with a (9) unrestricted in sign or magnitude 
p2|,p[. The fact that this behavior is not a limit of the generic AVTD solution (Q) illustrates further the qualitative 
distinction between the polarized case and Q ^ but arbitrarily small. 



III. NUMERICAL METHODS 



The numerical method is discussed in detail elsewhere [ |L However, we shall summarize it here. The work reported 
here was performed by using a symplectic PDE solver [ pl| , 22j . Consider a system with one degree of freedom described 
by q(t) and its canonically conjugate momentum p(t) with a Hamiltonian 



H 



£- + V(q)=H K +H v . 
2m 



(10) 



Note that the subhamiltonians Hk and Hy separately yield equations of motion which are exactly solvable no matter 
the form of V. Variation of Hk yields q = p/m, p — with solution 
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p(t + At)=p(t) , q(t + At) = q(t) + ^-At. (11) 

m 



Variation of Hy yields q — 0, p — —dV/dq with solution 

q(t + At) = q(t) , p(t + At) = p(t) 



dV 
dq 



At. (12) 



Note that the absence of momenta in Hy makes (|lj) exact for any V(q). One can then demonstrate that to evolve 
from t to t + At an evolution operator U( 2 )(At) can be constructed from the evolution sub-operators UK(At) and 
Uy(At) obtained from (0) and (|l2|). One can show that Q 

U {2) (At) = U K {At/2)U v {At)U K {At/2) (13) 

reproduces the true evolution operator through order (At) 2 . Suzuki has developed a prescription to represent the full 
evolution operator to arbitrary order p3| . For example 

%) (At) = U m (sAt) U {2) [(1 - 2s) At] U m (a At) (14) 

where s = 1/(2 — 2 1 / 3 ). The advantage of Suzuki's approach is that one only needs to construct U(2) explicitly. bt(2n) 
is then constructed from appropriate combinations of 

The generalization of this method to N degrees of freedom and to fields is straightforward. In the latter case, 
— * VW(%it)] so that dV/dq becomes the functional derivative SV/Sq. On the computational spatial lattice, 
the derivatives that are obtained in the expression for the functional derivative must be represented in differenced 
form. We note that, to preserve nth order accuracy in time, nth order accurate spatial differencing is required. Some 
discussion of this has been given elsewhere || . 

Convergence studies have been performed to test the algorithm. We have already shown || that at any fixed value 
of r, there exists a spatial resolution at which every spiky feature is resolved. Any finer resolution will yield identical 
results. However, as we shall see, generic spiky features grow narrower as the singularity is approached. This means 
(see [|J ) that a spatial resolution which was adequate at some To will fail to resolve all features at some t\ > tq. This 
type of "resolution dependence" is well understood to be a consequence of the physical behavior of the model. If the 
physical system has features which become arbitrarily narrow as r — > oo, no resolution will be adequate everywhere 
for all time. 



IV. GOWDY PHENOMENOLOGY 



Here we report the results of a numerical study of the evolution of generic Gowdy models toward the singularity 
(t — * oo). Rather than perform simulations for hundreds of Gowdy models, we consider a single class of models and 
argue that the observed behavior is generic. For the most part, we shall use the initial data P = 0, up — uocos#, 
Q = cos 9, and ttq — 0. This model is actually generic for the following reasons: The cos dependence is the smoothest 
nontrivial possibility. Since we shall discuss the growth of small scale spatial structure, we expect the cleanest effects to 
arise from the smoothest initial data. A more complicated initial state will not yield any qualitatively new phenomena. 
For example, choosing cos n9 instead yields the same solution repeated n times on the grid thus yielding the same 
result with poorer resolution. In addition, the amplitude of Q is irrelevant since the Hamiltonian (^J) is invariant 
under Q — > pQ, P — > P — In/? for any constant p. This also means that any unpolarized model is qualitatively 
different from a polarized (Q = 0) one no matter how small Q. Traveling waves (subject to J d9 pg = where pg 
is the total 0- momentum) yield qualitatively similar behavior 0). We are interested in the spatial structure that 
can develop from smooth initial data. For this reason, we have used cos# spatial dependence. Starting with a more 
complicated wave form will yield an eventually more complicated spatial structure but nothing qualitatively different 
will appear. Similarly, since rescaling Q does not change the solution qualitatively, the amplitude of Q is kept fixed. 
In this discussion, we shall consider only P and Q and construct A later if desired. 

From the chosen class of initial data, with high spatial resolution, one obtains after some time (r « 12) the profiles 
for P and Q shown in Fig. [l| The numbered features are typical for any generic solution. The peaks labeled #1, #2, 
and #3 arise from the same features of the nonlinear interactions. On a finer spatial scale in Fig. we notice that (a) 
the peak in P is associated with an extremum of Q (which may be a maximum or minimum) and (b) \ttq\ becomes 
large with ttq < (> 0) for a maximum (minimum) in Q. The apparent discontinuity in Q (labeled #4 in Fig. |l|) is 
shown on a finer scale in Fig. |^. We see that this feature occurs as ttq goes through zero where P < 0. In Fig. QTthe 
same features #1 and #4 are shown at r « 18. We see that the "spikes" have narrowed in 9. 
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The evolution of P and Q is shown in Fig. || (see also ||) where it is clearly seen that ever shorter wavelength 
modes develop until the point at which the AVTD behavior sets in. Fig. [| shows an interesting self-similarity in the 
early evolution. The parameter vo in the initial data is varied in a sequence of simulations. In each case, the number 
of peaks in P is counted after each time step. We find 

= CLN («o - vo,n) (15) 



where tjv is the time at which TV peaks are first present in the waveform and ajv is a constant. From (|15|), vo denotes 
the value of vq such that the Nth peak appears only at r = oo. The scaling shown in Fig. ^|is for N = 5. Similar (but 
less accurate) straight line fits are found for N equal 3 and 7. (The wave equations are even in 9 as are the initial 
data so that P and Q remain even functions of 9. Thus, except for a central peak, other peaks will form two at a 
time.) 

Thus, one may summarize that spatial structure forms more rapidly and evolves a more complicated form, the 
greater the initial "energy" in the waves as denoted by, e.g., v 2 from (||). Eventually, the waveforms of P and Q 
"freeze" to the AVTD behavior of P increasing linearly and Q constant in r. As has been emphasized by Hern 



and Stuart (HS) |25|, the duration in r of the observed non-AVTD regime depends on the spatial resolution of the 
simulation. 

In Fig. the maximum value over the spatial grid, v max , of the AVTD parameter v is shown as a function of time 
t for two simulations with different spatial resolutions. The steps in v max occur when the location of the spatial 
point with the maximum value of v changes. Note that, for the lower resolution simulation, v max eventually falls 
below unity in support of the GM conjecture jlq ]. We also see that the higher resolution simulation starts to diverge 
from the lower resolution one, signaling the presence of narrow features that are unresolved by the coarser resolution. 
(The former was not run longer primarily because the Courant condition requires small time steps for high resolution 
studies.) This suggests (as confirmed by HS) that the higher spatial resolution simulation will take longer before 
it becomes AVTD everywhere. This is shown explicitly in Fig. [8| where irp > 1 vs r is plotted for three different 
spatial resolutions. Note that for two of the spikes, P, T is larger for the finest resolution. In Fig. |^, the difference 
between P/t and v is plotted vs. r. Since P — > vt as r — ► oo, we see the influence of the next (constant in r) term 
P a + ln[(l + cosVO/2], in the AVTD solution for P (0). 



V. STRUCTURE IN THE GOWDY MODEL 



To understand the Gowdy phenomenology of Section III, we must explore the structure of the wave equations for 
P and Q. Let us consider (0) first. The equation for Q may be rewritten in first order form as 

Q,r = e^ 2P TT Q , 

n Q , T =e- 2T (e 2P Q, e ),o . (16) 

This immediately provides the explantion for both the large central peak which quickly forms in Q and for the apparent 
discontinuity in Q seen in Figs. [I] and [j[ From (|l6|), Q will grow rapidly where P is negative. Since, early in the 
evolution, P » vqt cos 9, one expects a large peak (as is observed in Figs, [j] and ||) where cos 9 < 0. We shall see later 
that P is driven to positive values eventually stopping the growth of Q. However, we shall also see that if ttq sa 0, P 
can remain negative longer. Say P w — Pq for Pq > at a value of 9 = 9% where ttq w ttq (9 — 9% ) . Thus we find that 

Q, T ttw° Q (9-9 1 )e 2P ° (17) 

so that \Q, T | is large and Q increases (decreases) exponentially depending (say for ttq > 0) on whether 9 > (<)#i- 
Note that we assume ttq, t rs 0. However, (|l6|b) implies this for P large and negative. 

More interesting is the equation for P. Consider two approximate forms of (0) depending on which nonlinear term 
dominates. Since the nonlinear terms depend exponentially on P, they will dominate the linear spatial derivative 
term when the exponential is large. 

For example, if P is large and negative, 

P, TT -tt 2 q e- 2P w (18) 

which has the first integral 

f, T +KQe =k v (19) 
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This allows one to define an effective potential 



(20) 



which is important when P, T < and (especially if) P < 0. This allows immediate explanation of the following 
phenomena: (a) In the AVTD regime, P > and its asymptotic time derivative v > since if P, r < 0, a bounce off 
Vi will occur at which P, T — > — P, T so that P will eventually reach positive values, (b) However, if ttq w 0, V% « so 
that both P and P )T can remain negative for a long time. Non-generic behavior occurs at the isolated point 6 = d\ 
where V% vanishes since P and P, T will never be driven to positive values there. 

On the other hand, if (a) P is large and positive and (b) P, T > 1, the wave equation for P is given by 



P*1TT 



,2(P-r) 



Q,2wQ 



(21) 



(22) 



which is easy to integrate (for Z = P — r) as 

Thus if P, T > 1, the dynamics of P is dominated by an effective potential 

V 2 =Q, 2 e e^ p - T \ (23) 

Again, from (fl6|si), our implicit assumption that Q, T ~ is consistent. One sees immediately, that if P, T > 1 a bounce 
off Vi will occur with the result that Z, T — > — Z, T or 



P, r ^-(P, r -2). 



(24) 



Thus the existence of V2 provides a mechanism to drive P, T and thus asymptotically v to values below unity. In the 
actual simulation, the asymptotic regions for the scattering are not reached so that ( p4| ) cannot be used to compute 
the numerical value of the change in P, T with r. Again one sees that non-generic behavior can arise, in this case, at 
9 = O2 where Q,g = 0. At 6*2 the mechanism to drive v below unity is absent so that larger values of v are allowed. 
If V2 is non-vanishing but small, it will take a long time for the bounce to occur so that the AVTD limit will take a 
long time to arise. This also explains the peaks in P seen in Fig. [j] Where Q,g w 0, Vz is small so that P can increase 
without hindrance with a large positive value of P, T . 

For a more quantitative understanding of the small scale spatial structure, note that both equations (18) and (21) 
have closed form solutions. Equation (18) holds in the AVTD regime, where the solution is given in equation (7). 
The exceptional point 9\ is the point where cqsi]) = — 1. Near this point we have cos ip » -1 + (ip'y(e - 6*1)72 where 
ip' = ?/>, 0(6*1). It then follows that for large r 



P w P Q - vt + In 

Q~Qo - 



1 + (V') V OT (0-0i)74 

e~ PQ ip' (6-61) 
2e -2„r + (^,') 2 (6/-6»i) 2 /2 



(25) 



Thus there is a spiky feature in both P and Q at 9 = 9% as is shown in Fig. |l^. Furthermore, these features 
steepen exponentially with time. Because of this steepening, the spatial derivatives of P and Q become large at late 
times. However, the AVTD equations are formed by dropping spatial derivatives from the exact equations (2) and 
(3). Therefore, one might worry that the terms that have been neglected are not actually negligible near the peaks. 
However, using equation (||) one can show that near 9 = 9\ the neglected terms go like e 2 '- 1 ' -1 ) 7 ' and are therefore 
exponentially damped for v < 1. 
The solution of equation (21) is 



P 



(26) 



Pq + t — In [cosh hit — cos <fr sinh wt] 
Here, Pn, w > 0, and 4> are functions of 9, and = e~ p ° w sin0. Note that P — ► (1 — w)t for large r where cos cf> ^ 1. 
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The exceptional point 6*2 is the point where cos 0=1. Near this point we have cos < 
(j)' = ipfi{6\). It then follows that for large r 



1 



P~ P + (l + tfl)T 



In 



b'f e 2wT 1 



? 2 ) 2 /4 



6 2 ) / 2 where 



(27) 



G 



Thus the feature in P narrows exponentially with time as can be seen in Fig. [l0[ The spatial derivatives of P become 
very large. However, the neglected terms in the equations go like e 2 *- 10-1 -* 1 ". Since w — 1 < for P T > asymptotically, 
it follows that these terms become negligible after the last bounce. 

The differences in Figs. [z]and || seen by increasing the spatial resolution are also understood in terms of our analysis 
of structure growth. At higher spatial resolution it is more likely that a spatial grid point of the numerical simulation 
will be very close to a non-generic point where Q,g — 0. Thus the non-generic behavior will become more apparent and 
the time to reach the AVTD limit characterized by v < 1 everywhere will be longer. Clearly, the observed resolution 
dependence may be interpreted to be evidence for the existence of non-generic points. The apparent discontinuity 
in Q seen in Figs. |l| and || is evidence for the non-generic point in 9 where ttq — 0. Its effect may be seen as a 
resolution dependence of the shape of this feature. It becomes steeper and narrower with increasing spatial resolution 
for sufficiently large t (see Fig. |ll|) or for increasing r with fixed spatial resolution (see Fig. |^b). As in Fig. ^, we 
note that the spikes in Fig. || where v > 1 are associated with extrema in Q. This is shown in Fig. [l2| Note that the 
alignment is not perfect between the peak in P, T and the extremum in Q where the spike is well resolved (see also 
Fig. ||a). This is understood as evidence that the AVTD regime has not yet been reached. In particular, the site of 
the extremum may be distorted by a nearby larger feature in Q. Running the simulation twice as long (Fig. ^) shows 
much closer alignment and is evidence that the predicted non-generic behavior will eventually dominate the spikes. 

The effective potentials V\ and V 2 are displayed in Fig. [L4|. We see that the approach to the AVTD regime is 
characterized by repeated alternate interactions of P with the two potentials. This is displayed explicitly for a typical 
value of 9 in Fig. [15] where P, V\, and V% are shown as functions of r. After each interaction with V2, the quantity 
P, T —1 changes sign. If this yields < P, T < 1, no further bounces occur. If P lT < 0, there will be an interaction 
with V\. Eventually, \P, T | falls below unity so that V 2 permanently disappears. P then continues to increase for all 
subsequent r with constant slope < v < 1 as is characteristic of the AVTD limit. (Note that V\ will also decrease 
exponentially as P increases even though it is part of the AVTD solution.) Examination of ([l6]) shows that, when 
\P, T I falls below unity and r — > 00, 7Tq, t — > as required by the AVTD limit. Once ttq becomes constant, as P 
increases, Q, T — > 0, again as required. The development of Q with time is shown in Fig. [lfi] for the same value of 
9 as in the previous figure. Here we see clearly how Q grows fastest when V\ dominates. This becomes even more 
dramatic in Fig. ^ which is the same graph for 9 ~ #1, the site of the apparent discontinuity in Q. 

It is important to emphasize that the magnitude and time evolution of V\ and V% as well as the initial values of the 
wave amplitude time derivatives are different at different values of 9. This explains the development of complicated 
small scale spatial structure. From ( [l9| ) or (|22|), one computes the time to the first bounce off V\ or V2 as 

r-H^/Kj) dx 
At= / (28) 

J ° Kjy/l-Vj/K* 

where J = 1 or 2 and Vj=VjsXt — 0. Here \i — ~P an d Xi = % ■ Evaluation of (|l^) and (p2|) at r = yields an 
approximate inverse scaling with vq seen in Fig. [| The central peak appears first and is easy to understand. Since, 
initially, P, T = focosf?, the highest negative velocity occurs at 9 — ir causing the bounce off V\ to occur there first. 
P then begins to increase at 9 = tt both earlier and faster than at neighboring values of 9. More complicated is the 
next pair of peaks which occur near 9 = 7r/4, 7tt/A. Initially, the profile of P vs 9 steepens. However, bounces off V2 
will begin to occur. This causes the tangent to P vs 9 to evolve from negative (positive) to positive (negative) slope 
for < 9 < it/ 2 (3ir/2 < 9 < 2ir) eventually creating a pair of peaks. The higher the original value of vo, the greater 
will be the number of bounces before V2 disappears for good. Ever finer small scale structures are created at 9 values 
where many bounces occur since P will increase for some spatial regions and decrease for others in a manner that 
changes with time. As v falls below unity, at a given 9, the bounces, and thus structure formation, will cease there 
but structure formation will continue in the ever smaller spatial regions where v continues to exceed unity. 



VI. DISCUSSION AND CONCLUSIONS 

We see that both potentials V\ and V2 drive the cosmology toward the AVTD regime characterized by (0) as t — * 00 
with < P, T — * v < 1. This is because Vi appears and influences the dynamics only when and where P, T < while 
V2 does so only for P, T > 1. In the AVTD regime, with < v < 1, both corresponding terms in the wave equation 
for P are exponentially small. We note that the wave equations (||) and (||) provide no other mechanism to drive P, T 
into the required range if it is not there initially. Recall that the polarized Gowdy model can have any value for P, T . 
Thus we argue that at isolated points (a set of measure zero in 9) where cither ttq or Q,g is precisely zero, P, T could 
remain outside the allowed range. In the former case (say at 9 = 61), the consistency equations (^|) are satisfied so 
that one expects the solution at 6\ to be like an AVTD polarized solution. However, in the latter case (say at 9 = 9 2 ), 
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the consistency conditions require Q,gg = as well as Q,g = — something one does not expect generically. If Q,gg = 
at 02, then the solution is of the AVTD large v type there. If Q : gg^ 0, the behavior cannot be analyzed either 
numerically or by these analytic methods. We may therefore conclude that the simulations provide strong support for 
the conjecture that the Gowdy models have a singularity which is AVTD except, perhaps, at a set of measure zero. 

The existence of this set of measure zero of "non-generic" spatial points can be inferred from the simulations since 
near them the potentials are very flat. Thus, the closer one is to e.g. Q\ or 62 defined as above, the longer it will take 
to drive P, T to the range [0, 1) there. This manifests itself as a resolution dependence of the simulations of a very 
specific type. At any value of r (say f), there is a spatial resolution sufficient to resolve all spiky features. (This was 
demonstrated in ||.) But the "center" of each spiky feature is the site of a non-generic point. As one moves away 
from this point in either direction in 9, the coefficient of the potential will be larger. Thus the interaction with the 
potential will occur first at regions away from the center of the feature and move closer to the center as r increases. 
Since the effect of the interaction is to change the sign of either P, T or P, T — 1, the interactions will tend to eliminate 
the outer part of the spiky feature. Thus the spiky feature gets narrower in time until for some r > r , the original 
resolution is no longer adequate. Although, in a peak near Q,g = 0, P, T should be > 1, inadequate resolution may act 
as an effective averaging to give a measured smaller value of P, r . Thus the peak in v appears to decrease in amplitude 
only because, at most 9 values, the interaction with V2 has already occurred. We note that this effective averaging 
will give < v < 1 everywhere for some rp for any fixed spatial resolution. As the spatial resolution is increased, 
tf will also increase since accurate resolution of spiky features can be maintained farther into the simulation and the 
effective averaging will take longer to come into play. Thus the apparent resolution dependence is consistent with the 



behavior of the potentials V\ and V2 ex- 
Several examples are shown in Figs 



cpected near the non-generic points. 
18|-|2l]. In Fig. |lj, a spiky feature in P, T is shown for three different spatial 
resolutions at r = 24.79. While the coarsest resolution shows v < 1, it also clearly fails to resolve the spiky feature. 
The highest resolution both resolves the spiky feature and shows that it has a core where v > 1 as expected. A similar 
spiky feature in P, T is shown in Fig. |l9| for three different times for the same spatial resolution. One expects that the 
clearly resolved spike observed early in the simulation will narrow as r increases remaining always > 1 at the center of 
the peak. This is not observed (rather one sees the constant P, T behavior expected in the AVTD limit). However, one 
again sees clearly that the spatial resolution which was adequate at r = 12.4 has become inadequate at later times 
and that the properties of the spike are not correctly represented. Fig. ^ shows that the apparent discontinuity in 
Q becomes narrower and steeper with time and eventually is not adequately resolved at a given spatial resolution. 
One expects P to be increasingly negative at the center of this feature since P, T < is maintained at that non-generic 
point. At fixed spatial resolution, we see in Fig. |2l] that the negative region of P both narrows and deepens as r 
increases from 12.4 to 18.6 as expected. However, at r = 24.8, the expected behavior is not seen since the depth of the 
negative region fails to increase. However, we also see that the spatial resolution has become inadequate. Presumably, 
adequate spatial resolution would show the narrowness and depth of the negative region suggested by the curvature 
of P vs 9 in the region just outside the core. 

In conclusion, the behavior of generic Gowdy cosmologies on T 3 x R can be completely understood in terms of 
the nonlinear interactions between the two polarizations of the gravitational waves. These act as effective potentials 
which drive the system to the prediced AVTD behavior and then cease to play a role. At the set of measure zero 
where the potentials are absent, P, T may lie outside its allowed AVTD range of [0, 1) for all r if it does so at any 
time since the mechanism to correct its value is absent. The existence of this set of measure zero is inferred from the 
details of the shape and time dependence of the spiky features and from the inability to resolve them with a given 
spatial resolution beyond a certain time. 



ACKNOWLEDGEMENTS 



The authors would like to thank Vincent Moncrief and Boro Grubisic for useful discussions. BKB would like to 
thank the Institute for Geophysics and Planetary Physics at Lawrence Livermore National Laboratory and the Albert 
Einstein Institute at Potsdam for hospitality. This work was supported in part by National Science Foundation 
Grants PHY9507313 and PHY9722039 to Oakland University. Computations were performed at the National Center 
for Supercomputing Applications (University of Illinois). 



[1] V.A. Belinskii, E.M. Lifshitz, I.M. Khalatnikov, Sov. Phys. Usp. 13, 745-765 (1971); Adv. Phys. 19, 525-573 (1970). 



8 



D. Eardley, E. Liang, R. Sachs, J. Math. Phys. 13, 99 (1972). 

J. Isenberg, V. Moncrief, Ann. Phys. (N.Y.) 199, 84-122 (1990). 
C.W. Misner, Phys. Rev. Lett. 22, 1071-1074 (1969). 

M.P. Ryan, Jr., L.C. Shepley, Homogeneous Relativistic Cosmologies (Princeton University, Princeton, 1975). 

E. Kasner, Am. J. Math 43, 130 (1921). 

J.D. Barrow, F. Tipler, Phys. Rep. 56 372 (1979). 

B.K. Berger, V. Moncrief, Phys. Rev. D 48, 4676 (1993). 

B.K. Berger, in Proceedings of the 14th Internatial Conference on General Relativity and Gravitation, edited by M. Fran- 
caviglia, G. Longhi, L. Lusanna, E. Sorace (World Scientific, Singapore, 1997). 

B.K. Berger, D. Garfinkle, V. M oncrief, in Proceedings of the Workshop on Black Hole Interiors, edited by L. Burko, A. 



Ori (Haifa, 1997), gr-qc/9709073 



R.H. Gowdy, Phys. Rev. Lett. 27, 826-829 (1971). 

B.K. Berger, Ann. Phys. (N.Y.) 83, 458 (1974). 

V. Moncrief, Ann. Phys. (N.Y.) 132, 87-107 (1981). 

A. Einstein, N.J. Rosen, J. Franklin Inst. 223, 43 (1937). 

B. Grubisic, V. Moncrief, Phys. Rev. D 47, 2371 (1993). 

P.T. Chrusciel, J. Isenberg, V. Moncrief, Class. Quant. Grav. 7, 1671 (1990). 

Ch. Charach, S. Malin, Phys. Rev. D 21, 3284 (1980). 

M. Carmeli, Ch. Charach, S. Malin, Phys. Rep. 76, 79-156 (1981). 

M.P. Ryan, Jr., "Plane-Symmetric Cosmologies — An Overview," in SILARG, pp. 96-121. 

B.K. Berger, V. Moncrief, "Numerical Evidence for a Velocity Dominated Singularity in U(l) Symmetric Cosmologies," 
unpublished. 

J. A. Fleck, J. R. Morris, M.D. Feit, Appl. Phys. 10, 129-160 (1976). 

V. Moncrief, Phys. Rev. D 28, 2485-2490 (1983). 

M. Suzuki, Phys. Lett. A146, 319 (1990). 

CM. Swift, Masters Thesis, Oakland University (1993). 

S.J. Hern, J.M. Stewart, "The Gowdy T 3 Cosmologies Revisited," |gr-qc/97080"3^. (However, see also ^r-qc/97080"5^.) 



120 



80 
40 

Q 





-40 

1 q 2 3 

FIG. 1. P (solid line) and Q (dashed line) vs 9 at r = 12.4 for the standard initial data set, P = 0, np = Vo cos#, Q — cos#, 
7tq = 0, with vo = 5 for < 9 < n. The simulation was run with 20000 spatial grid points in the interval < 9 < 2-n. The 
numbers on the graph refer to the most interesting features. Peaks #1, #2, and #3 in P are essentially the same in that they 
occur at extrema (in 9) of Q. (Q has a maximum at the locations of peaks #1 and #3 and a minimum at #2.) Structure #4 is 
qualitatively different in that there appears to be a discontinuity in Q. This occurs where P<0 (P = 0is shown as a dashed 
line) and ttq ~ 0. 
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FIG. 2. Details of feature #1 in Figure (a) P (solid line) and Q (dashed line) vs 9 near the peak. Note the scales in 
and 9. (b) P (solid line) and ttq (dashed line) vs 9 near the peak. 
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FIG. 3. Details of feature #4 in Figure (a) P (circles) and ttq (solid line) vs 9. The horizontal dashed line indicates ttq = 
while the vertical dashed line marks the 9 value at which the ttq goes through zero. Note the offset between the minimum in 
P and the zero of ttq. (b) ttq (solid line) and Q (dot-dashed line) vs 9 near the feature. 
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FIG. 4. Evidence for the narrowing of spiky features with increasing r. (a) P as shown in Figure ^| is plotted at two different 
times, (b) Q as shown in Figure ^ is plotted at two different times. The amplitudes of P and Q have been rescaled for clarity. 




FIG. 5. P, Q, X, and U = e~ 2r (Pj + e 2P Q,f) are shown in the 0-t plane. The arrow is in the direction of increasing r. 
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FIG. 6. Scaling in the Gowdy Model. Plot of 1/ts, the inverse of the time at which the 5th peak appears in P vs vo- Two 
cases are shown for initial data P — 0, Q = cos8: (1) irp = 1/V2 vocosfl, ttq = l/y/2 vocosfl is indicated by filled circles; (2) 
ttp = vo cosd, 7vq — is indicated by open squares. The solid line is a linear best fit. 
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FIG. 7. v m ax vs t for the Gowdy model with vo = 5. The solid curve is from a simulation with 3200 spatial grid points 
while the open circles represent one with 20 000 spatial grid points. The horizontal line is vo = 1. The difference near the end 
of the higher resolution simulation is due to the influence of the non-generic point where Q,e = 0. 
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FIG. 8. P, T > 1 vs e at r 
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24.79 with standard initial data with vo = 10 for three spatial resolutions. For the coarser 



resolutions, all simulation spatial points are shown. 
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FIG. 9. The same as Figure ^ with the maximum value of P/r — v over the spatial grid also plotted (dashed line) vs r. 
observed behavior is due to the form P — > vt + ln[(l + cos^)/2] in the AVTD limit asr-» oo. 
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FIG. 10. Sketches of the analytic approximations to the peaks in P and Q. The larger r value is indicated by the broken 

line, (a) P from (25a). (b) Q from (25b). (c) P from (27). 
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FIG. 11. Q vs 6 for r = 24.79 with vo = 10 at the site of the apparent discontinuity for three spatial resolutions. At this 
value of r, the feature is not resolved even at the finest resolution used. 
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FIG. 12. Details of P (solid line, left axis scale) vs Q (dot-dashed line, right axis scale) for the spikes in Figure H (a) Note 
that the spike is completely resolved and that there is an offset of the extrema. The inset shows that this peak is part of a 
larger feature in Q. (b) Again the offset can be related to a larger feature in Q. (c) The extrema are much more closely aligned, 
(d) The features are not completely resolved. The extrema are aligned to within the resolution, (e) Again the features are not 
well resolved. The change in the value of Q is the smallest that can be seen in single precision. The simulations were done in 
double precision but the data was reported in single precision. 
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FIG. 13. Spiky features (a) and (b) of P (solid lines, left axis scale) with the associated Q (dot-dashed lines, right axis scale) 
in Figure ^ at r = 49.5. Part of the previous figure (at r = 24.8) is shown with finer lines. Note the improvement in alignment 
and the narrowing of the spike in P. 
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FIG. 14. V\ vs P and Vi vs P — r shown as solid lines. The horizontal dashed arrows denote the constants n\ and K2- 
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FIG. 16. Q (solid curve), V\ (dash-dotted curve), and V2 (dashed curve) vs r at the same value of 9 as in Figure FT 




FIG. 17. Q (solid curve), Vi (dash-dotted curve), and — P (dashed curve) vs r at the site of feature #4 in Figure 
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FIG. 18. The same spiky feature in P, T vs 8 is shown at the same r for resolutions of 1024 (dash-dotted line), 4096 (dashed 
line), and 40960 (solid line) spatial grid points. The horizontal line marks v = 1. 
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FIG. 19. For a fixed resolution of 20224 spatial grid points, a spiky feature in P, T is shown vs 9 for three different values of 
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FIG. 20. The apparent discontinuity in Q vs 9 at three values of r with 20224 spatial grid points. 
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FIG. 21. Part of the region where P < vs 9 at three values of r for 20224 spatial grid points. 
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